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Uniqueness problems in the elliptic sector of constrained formulations of Einstein equations have 
a dramatic effect on the physical validity of some numerical solutions, for instance, when calculating 
the spacetime of very compact stars or nascent black holes. The fully constrained formulation (FCF) 
proposed by Bonazzola, Gourgoulhon, Grandclement, and Novak is one of these formulations. It 
contains, as a particular case, the approximation of the conformal flatness condition (CFC) which, 
in the last ten years, has been used in many astrophysical applications. The elliptic part of the FCF 
basically shares the same differential operators as the elliptic equations in CFC scheme. We present 
here a reformulation of the elliptic sector of CFC that has the fundamental property of overcoming 
the local uniqueness problems. The correct behavior of our new formulation is confirmed by means 
of a battery of numerical simulations. Finally, we extend these ideas to FCF, complementing the 
mathematical analysis carried out in previous studies. 



PACS numbers: 04.20.Ex, 04.25.Nx, 04.25.D-, 97.60.-s 

I. INTRODUCTION 

In recent years we have seen the successful application 
of numerical codes to accurately calculate the spacetimcs 
of compact astrophysical objects like collapsing stellar 
cores, (proto)neutron stars, and black holes. Most of 
these codes are based on the 3 + 1 formalism of general 
relativity (see, e.g., P, 0, [l] for reviews). They typically 
fall into two classes. One approach relies on the free 
evolution of the 3 + 1 Einstein equations, recast in order 
to cure long-term stability problems. Here the constraint 
equations are only solved initially, and closely monitored 
at each time step to control the accuracy of the numerical 
solution. 

Alternatively, formulations based on a constrained evo- 
lution, where the constraints are solved in parallel with 
evolution equations, have proven to be successful as well. 
Such approaches exhibit the advantage that the solution 
cannot violate the constraints by definition (within the 
accuracy of the numerical scheme). In particular, the 
conformally flat approximation [1, Q (hereafter CFC) of 
the full Einstein equations, which constitutes a fully con- 
strained formulation, has been shown to yield long-term 
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stable evolutions of such astrophysical scenarios (see, e.g., 
[1, S B 0)- However, apart from computational chal- 
lenges, arising from the need to frequently solve the el- 
liptic constraint equations, constrained formulations suf- 
fer from mathematical nonuniqueness problems when the 
configuration becomes too compact. In the case of the 
collapse of a stellar core or a (proto)neutron star to a 
black hole, such a situation is encountered already be- 
fore the apparent horizon forms. This issue has, in the 
past, been prohibitive to successfully applying such for- 
mulations in numerical simulations of a wide range of 
astrophysical problems. 

The nonuniqueness of solutions stems from the non- 
linearity of the constraint equations and has been stud- 
ied within the so-called extended conformal thin sand- 
wich (XCTS) [l^, [111, approach to the initial data 
problem in general relativity. In Ref. (isj a parabolic 
branching was numerically found in the solutions to the 
XCTS equations for perturbations of Minkowski space- 
time, providing the first evidence of nonuniqueness in 
this elliptic system. First analytical studies have been 
carried out in [T^ . Il5| . finding support for the gener- 
icity of this nonuniqueness behavior. More specifically, 
the XCTS elliptic system is formed by the Einstein con- 
straint equations in a conformal thin sandwich ( CTS) 
decomposition [loj supplemented with an additional el- 
liptic equation for the lapse function, which follows from 
the maximal slicing condition. Although no general re- 
sults on existence and uniqueness for the XCTS system 
are available (in contrast to the CTS case and similar 
elliptic systems encompassing only the constraints; see. 
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e.g., [ll,[i3,[ni,[ll,[i3)' the analysis in ^ strongly sug- 
gests the presence of a wrong sign in a certain term of the 
lapse equation as the culprit for the loss of uniqueness, es- 
sentially because it spoils the application of a maximum 
principle to guarantee uniqueness. Moreover, in these cir- 
cumstances (namely, the existence of a nontrivial kernel 
for the XCTS elliptic operator) it is shown in [l^ that 
the parabolic behavior found in is indeed generic. 

Certain constrained evolution formalisms which incor- 
porate elliptic gauges in their schemes contain elliptic 
subsystems which share essential points with the XCTS 
equations. Nonuniqueness in the elliptic subsystem is 
certainly an issue for the well-posedness of the whole 
elliptic-hyperbolic evolution system. In numerical im- 
plementations this can depend on the employed numer- 
ical scheme, in particular, on its capability to remain 
close to one of the solutions, at least as long as the so- 
lution stays sufficiently far from the branching point. In 
fact, constrained or partially constrained evolutions have 
shown to be robust in a variety of contexts (see, e.g., 
the references in [31 and Sec. 5.2.2 of [l^). However, 
the problems described above have also emerged, for in- 
stance, in the axisymmetric case in [20,[2l| (see also (2^). 
The analysis in concludes that the reason behind the 
failures in these axisymmetric formulations is in fact re- 
lated to the presence of wrong signs or, more precisely, to 
the indefinite character of certain non-linear Helmholtz- 
like equations present in the scheme (see [11] for details 
and also for a parallel numerical discussion in terms of a 
class of relaxation methods for the convergence of the 
elliptic solvers). Regarding the full three-dimensional 
case, fully constrained formalisms have been presented 
in [ii, [11 [ii]. While the work in [13, ^ includes an 
elliptic subsystem closely related to the XCTS equations 
and therefore suffers potentially from these nonuniquc- 
ness problems, the uniqueness properties of the scheme 
of [2^ must yet be studied. In both cases, the full nu- 
merical performance still has to be assessed. 

The goal of the present work is to discuss a scheme ad- 
dressing the nonuniqueness issues of XCTS-like elliptic 
systems in the full three-dimensional case, with astro- 
physical applications as our main motivation. Having 
the analysis of the fully constrained formalism (hereafter 
FCF) of [H, as our ultimate aim, we focus on an 
approximation in the spirit of the CFC approximation 
by Isenberg, Wilson, and Mathews [3, HE]. This method- 
ological choice is justified since the CFC scheme already 
contains the relevant elliptic system of FCF, but in a 
setting in which potential additional problematic issues 
related to the FCF hyperbolic part do not mix up with 
the specific problem we are addressing here. Therefore, 
we discuss in detail a modification of the CFC scheme (in 
the presence of matter) where maximum-principle lines 
of reasoning can be used to infer the uniqueness of the 
solutions. We investigate numerically the performance of 
the new CFC scheme and finally indicate the main lines 
for its generalization to the full Einstein FCF case. 

The article is organized as follows. In Sec. |lT] we re- 



view the FCF and CFC formalisms, and then discuss 
the limitations found in the numerical implementations 
of the latter. In Sec. IIIII we introduce the modification 
of the CFC scheme, with the aim of solving the unique- 
ness issues, and we present various numerical tests of the 
new scheme in Sec. IIVI In Sec. |V]the guidelines for the 
generalization to the FCF case are discussed, and conclu- 
sions are drawn in Sec. IVII In the Appendix we justify a 
further approximation assumed in Sec. IIIII which is con- 
sistent with the CFC setting. Throughout the paper we 
use the signature (—,-!-,+,+) for the spacetime metric, 
and units in which c = G = AIq = 1. Greek indices run 
from to 3, whereas Latin ones from 1 to 3 only. 

II. THE FULLY CONSTRAINED FORMALISM 
AND THE CONFORMAL FLATNESS 
CONDITION 

A. A brief review of the fully constrained 
formalism 

Given an asymptotically flat spacetime (A^jg^j/) we 
consider a 3 -I- 1 splitting by spacelike hypersurfaces Ej, 
denoting timelike unit normals to St by n^. The data 
on each spacelike hypersurface are given by the pair 
{'jijjK''^), where 7^^ = g^i, + n^riv is the Riemannian 
metric induced on Ej. We choose the convention K^,^ ~ 
— ^Cnjf_ii^ for the extrinsic curvature. With the lapse 
function N and the shift vector /3* , the Lorentzian metric 
gfj^i, can be expressed in coordinates (x^) as 

g^,^ dx^ dx" = -N^ dt^ + j,j (dx' + f3' dt) {dx^ + [3' dt). 

(1) 

On the other hand, we can write 

2NK^3 = dtf^ -f D^(3^ + D^I3\ (2) 

where Di is the Levi-Civita connection associated with 
7^^ and (?t7*-' represents the Lie derivative with respect 
to the evolution vector t^' {dtY = Nn^^ + (3^^. As 
in [2^ we introduce a time independent flat metric /y , 
which satisfies Ctfij = dtfij = and coincides with 7^ at 
spatial infinity. We define 7 := det7y and / := detfij. 
This fiducial metric permits the use of tensor quantities 
rather than tensor densities. The next step in the formu- 
lation of [2^ is the conformal decomposition of the 3 4-1 
fields. First, a representative 7^ in the conformal class 
of 7ij is chosen, so we can write 

7., = i^Hj , K^' = i^'^-^A^' + ^KY' , (3) 

where K = Kij and 7 := det7y , and C G IR- In 
Ref. [2^, the choice C = 4 was adopted, leading to the 
following expression of A'-' in terms of the lapse N and 
shift 

A^' = ^ (0^3^ + D^(3^ - lokP^f' + daA , (4) 
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Di being the Levi-Civita connection associated with . 
This is in the spirit of the decomposition employed in the 
(X)CTS approach to initial data. Regarding the choice 
of the representative of the conformal metric ja, a uni- 
modular condition 7 = / was adopted in [23| . so that 
ijj = (7//)^^^^- The deviation of the conformal metric 
from the flat fiducial metric is denoted by /i*-' , i.e. 

h'^:^f^^fi_ (5) 

Once the 3 + 1 conformal decomposition is performed, a 
choice of gauge is needed in order to properly reformulate 
the Einstein equations as partial differential equations. 



At/; 



A{NiP) 



A/3* + \v'VjP^ 



where A stands for the flat Laplacian (A := f^^ViVj), E, 
5* and S are, respectively, the energy density, momentum 
density, and trace of the stress tensor, all measured by the 
observer of 4- velocity n'^ (Eulerian observer) : in terms of 
the energy-momentum tensor T^j,, E T^^rT.^ri'^, S*' 
-7*T^,n^ and S := -f'^S^,, with S,, T^.7^»7"j- 
Furthermore, 

{Lpy^ := V'(3^ + V^p' - p'WkP", (11) 
A^- := ^7'' + T^.la ~ ^^a^,) • (12) 

Equation ([7]) follows from the Hamiltonian constraint, 
whereas Eq. ^ results from the momentum constraint 
together with the preservation of the Dirac gauge in time. 
Equation ([5]) corresponds to the preservation in time of 
the maximal slicing condition, dK/dt = 0. Note that 
expression ([TU]) for the Ricci scalar of the conformal met- 
ric does not involve any second order derivative of the 
metric; this property follows from Dirac gauge [2^. The 



The prescriptions in [23l | are maximal slicing and the so- 
called generalized Dirac gauge, 

K = 0, Vkf' = 0, (6) 

where 1?^ stands for the Levi-Civita connection associ- 
ated with the flat metric fij . The Einstein equations then 
become a coupled elliptic-hyperbolic system to be solved 
for the basic variables /i*^ , ip, N, and /3* [2^. 

Expressing the differential operators in terms of the 
connection of the flat metric, the elliptic part can be 
written as 



I 

resul ting elliptic subsystem coincides with the XCTS sys- 
tem jll| . except from the field chosen to solve the max- 
imal slicing equation: Eq. ([5]) above is to be solved for 
N^p, whereas in [ll| the conformal lapse N := Ntp~^ is 
employed instead. This directly affects the value (and, in 
particular, the sign) of the power of the conformal factor 
in the nonlinear terms of Eqs. ([7]) and ([5]). More gen- 
erally, one could define a generic rescaling of the lapse, 
N = Nip"-, such that the choice in [ll| corresponds to 
a = 6, whereas the choice in Eq. ([8]) above corresponds 
to a = — 1 (see [l^] for the general equations in the vac- 
uum case). An important remark is the absence of a 
choice of a such that the factors multiplying ip and N on 
the right hand side of the linearized versions of Eqs. ([7]) 
and ([S]) both present a positive sign. In the presence of 
matter, terms multiplying the energy density E also con- 
tribute to these sign difficulties, though in this case they 
can be fixed by an appropriate conformal rescaling of the 
energy density (see later). An additional concern in a 
generic evolution scenario is the sign of i?, which is also 
relevant in the linearized equations. Implications of this 
issue are discussed in Sec. IIIII 
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The Einstein equations in the form of the cUiptic equa- 
tions Q-dp and the hyperbohc equation for h^^ as given 
in Ref. [23| are to be solved together with the hydrody- 
namic equations, 



(13) 
(14) 



where is the Levi-Civita connection associated with 
the metric 5^^, p is the rest- mass (baryon mass) density, 
and is the 4- velocity of the fluid. 



B. The conformal flatness approximation 

If the hyperbolic part of the FCF system is not solved, 
but rather the condition ft,'-' = is imposed, the resulting 
three- metric jij is conformally flat, and the CFC approx- 
imation is recovered. Therefore, the FCF is a natural 
generalization of the CFC approximation. The latter has 
been used in many astrophysical applications, like the ro- 
tational collapse of cores of massive stars (6^, ,29|, [s^] or 
supermassive stars the phase-transition-induced col- 
lapse of rotating neutron stars to hybrid quark stars Q, 



equilibrium models of rotating neutron stars [31. 
well as for binary neutron star merger 0, HF 
The elliptic subsystem of the FCF, Eqs. ©-(El), reduces 
in CFC to 



E* 



16n 



E* + 2S* 
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(15) 



,(16) 
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N 



A/3'+-V'VjP^ = mTrN^j-^{S*y+ 2^P^°K'^Vj-^, (17) 

where the following rescaled matter quantities have been 
introduced, following York [H: 



E* 
S* 




(18) 
(19) 
(20) 



Equations psp and ()16p inherit the local nonuniqueness 
problems already present in the FCF equations. Al- 
though the sign problems specifically related to the en- 
ergy density terms are solved by the conformal rescaling 
of the components of the energy-momentum tensor and 
the CFC eliminates the R term, problems related to the 
KijK^^ term remain in the scalar CFC equations. This 
is apparent once the extrinsic curvature is expressed in 
terms of the lapse and the shift. 

Conformal rescaling of the hydrodynamical variables 
is not only relevant for local uniqueness issues. The hy- 
drodynamic equations (|13p and (fTH) can be formulated 
as a first-order hyperbolic system of conservation equa- 
tions for the quantities (D* , {S*)i, E*) [H, [si], where. 



similarly to Eqs. ((T8])-(l20l), D* := ^^D, D := Nu° p 
being the baryon mass density as measured by the Eu- 
lerian observer. We can thus consider E* and {S*)i as 
known variables in the computation of the CFC metric. 
Note that these quantities differ from E and Si by a 
factor -0^, and hence it is not possible to compute the 
nonstarred quantities before knowing the value of ip. If 
the energy-momentum tensor represents a fluid, then the 
source of Eq. (jl6[) cannot be explicitly expressed in terms 
of (D*, {S*)i,E*), the reason for that being the depen- 
dence of S* on the pressure P. The pressure can only be 
computed in terms of the "primitive" quantities, e.g., as 
a function P{p, e) of the rest-mass density and the spe- 
cific internal energy e. The primitive quantities are, in 
general, recovered from (D, Si, E) implicitly by means of 
an iteration algorithm. So far, two solutions of the prob- 
lem related to the fact that S* directly contains P have 
been used in numerical simulations performed with the 
CFC approximation. 

The first approach [2^ is to consider P, and hence also 
S* , as an implicit function of ip. Then Eqs. p^ -(fT7 |) can 
be solved as a coupled set of nonlinear equations using a 
fixed-point iteration algorithm. The convergence of the 
algorithm to the correct solution depends not only on the 
proximity of the initial seed metric to the solution, but 
also on the uniqueness of this solution. The latter point is 
extensively discussed in Sec. IIIII Furthermore, one prob- 
lem of this approach is the necessity of performing the re- 
covery of the primitive variables (which is numerically a 
time consuming procedure) to compute the pressure dur- 
ing each fixed point iteration. Because of the uniqueness 
problem, this approach can be only successfully applied 
in numerical simulations for, at most, moderately strong 
gravity (like stellar core collapse to a neutron star or the 
inspiral and initial merger phase of binary neutron stars) , 
but fails for more compact configurations like the collapse 
of a stellar core or a neutron star to a black hole. For 
such scenarios with very strong gravity, one finds conver- 
gence of the metric to a physically incorrect solution of 
the equations or even nonconvergencc of the algorithm. 

A second approach to the recovery algorithm problem 
is the attempt to calculate P independently of the CFC 
equations. This can be achieved by computing the con- 
formal factor by means of the evolution equation 



(21) 



The conformal factor tp' obtained in this way is analyt- 
ically identical to the tp from Eqs. (fT5|) -(fT7 |l . but here 
we use a different notation to keep track of the way it 
is computed. The value of tp' is solely used to evaluate 
P, and the coupled system of Eqs. P^ - (fT7|) is solved 
for determining ijj, N, and /3\ Although this approach 
allows one to avoid the problem of recovering the prim- 
itive variables at each iteration, it also suffers from the 
convergence problem, and the simulation of configura- 
tions with very strong gravity is still not feasible. Fur- 
thermore, new complications arc introduced by using two 
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difFcrcntly computed values, ip and ■0', of the same quan- 
tity. For some scenarios like the formation of a black 
hole from stellar collapse, the numerical values of these 
two quantities during the evolution of the system start 
to diverge significantly at some point. We find that this 
inconsistency cannot be avoided, since any attempt to 
artificially synchronize both values leads to numerical in- 
stabilities. 



since it involves the evaluation of the conformal factor. 
Nonunique solutions of "0, either due to the use of the 
nonconformally rescaled E or the quadratic extrinsic cur- 
vature term, spoil the convergence of the algorithm when 
density, and thus compactness, increases. We again em- 
phasize that a possible synchronization of and ip' does 
not solve the problem in general, since numerical insta- 
bilities eventually arise at sufficiently high compactness. 



III. THE NEW SCHEME IN THE 
CONFORMALLY FLAT CASE 

A. Uniqueness of the elliptic equations and 
convergence of elliptic solvers 

Well-posed elliptic partial differential systems admit 
non-unique solutions whenever the associated differential 
operator has a nontrivial kernel. When discussing suffi- 
cient conditions guaranteeing uniqueness, it is illustrative 
to first consider the case of a scalar elliptic equation. In 
particular, for the class of scalar elliptic equations for the 
function u of the form 

Au + huP ^ g, (22) 

where h and g are known functions independent of u, a 
maximum principle can be used to prove local uniqueness 
of the solutions as long as the sign of the exponent pis 
different from the sign of the proper function h [l], Is^. l38l. 

In the CFC case, we are not dealing with a single scalar 
elliptic equation, but rather with the coupled non-linear 
elliptic system (fT5|) - (fT7| . Therefore, assessing wheter or 
not the scalar equations (jlSp and present the good 
signs for the application of a maximum principle is an im- 
portant step for understanding the uniqueness properties 
of the whole system. However, as pointed out in the pre- 
vious section, the CFC equations for the conformal factor 
and the lapse possess the wrong signs in the quadratic ex- 
trinsic curvature terms (once everything is expressed in 
terms of the lapse and the shift). This problem can be 
fixed in Eq. (jl5|) by an appropriate rescaling of the lapse, 
N = Ntp^, but this strategy does not solve the problem 
for the lapse equation (cf. the discussion on the conformal 
lapse N in Sec. Ill A| ). Therefore, we cannot use the max- 
imum principle to infer local uniqueness of the solutions 
to the CFC equations. In these conditions of potential 
nonunique solutions, convergence to a undesirable solu- 
tion may happen. As mentioned in the introduction, this 
pathology has been illustrated using simple analytical ex- 
amples of scalar equations of the type (|22p in |14| , as well 
as in numerical implementations of the vacuum Einstein 
constraints in the XCTS approach and certain con- 
strained evolution formalisms (see, e.g., [111). 

In the context of the CFC approximation this sign 
issue has also appeared, in particular associated with 
the "recovery algorithm" problem discussed in Sec. Ill Bl 



B. Numerical examples 

The nonuniqueness of solutions has also been observed 
in FCF, as described in the following example. Let us 
consider a vacuum spacetime, with initial data formed by 
a Gaussian wave packet, as in [2^, but with much higher 
amphtude xo — 0.9 instead of xo = 10~^ in [11] (see the 
latter reference for notations). The integration technique 
and numerical settings are the same as in (23j . but con- 
trary to the results for small amplitudes obtained in that 
reference, the wave packet does not disperse to infinity 
and instead starts to collapse. Fig. [1] displays the time 
evolution of the central lapse A'c at r = and of the sys- 
tem's Arnowitt-Deser-Misner (ADM) mass Madm, which 
in the present conformal decomposition can be expressed 
as 

Madm ^ f (^'V- - dA' 

= -^j>V,^dA\ (23) 

where the integral is taken over a sphere of radius r = oo 
and the second equality follows from the use of Dirac 
gauge [Eq. (O]. 

The very sudden change at t ~ 0.4 in both the cen- 
tral lapse and the ADM mass, which is also present in, 
e.g., the central conformal factor Tpc, originates from the 
convergence of the elHptic system ([T])-® to another solu- 
tion with a different (unphysical) value of the ADM mass. 
The good conservation of A/adm and the smooth evolu- 
tion of Nc for t > 0.4 indicate that this other solution 
remains stable until i ~ 2, when high-frequency oscilla- 
tions appear. These oscillations may be due to the over- 
all inconsistency of the system, destabilizing the whole 
scheme. On the other hand, the time evolution of h^^ 
does not show any such type of behavior, and h^^ exhibits 
a continuous radial profile at all times. This is numerical 
evidence that, also for the full Einstein case (i.e. without 
approximation), the generalized elliptic equations suffer 
from a similar convergence problem as in the CFC case. 

The same subject is also exemplified when one tries 
to calculate the spacetime metric for an equilibrium neu- 
tron star model from the unstable branch using either 
Eqs. dll)-® in the FCF case or Eqs. (US])-!!!]) in the CFC 
approximation. Even for the simple setup of a polytrope 
with adiabatic index F = 2 in spherical symmetry, those 
metric equations yield - when converging at all - a grossly 
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FIG. 1: Time evolution of the central lapse A^c (top panel) 
and the ADM mass Madm (bottom panel) for a collapsing 
packet of gravitational waves, using the integration scheme 
proposed in [S^l . The unit of t is given by the initial width of 
the wave packet 



incorrect solution if the matter quantities [D,Si,E) in 
the source terms are held fixed. Both the metric compo- 
nents as well as the ADM mass can deviate from the phys- 
ical solution by a few tens of percent, even though that in- 
correct metric satisfies the asymptotic flatness condition. 
The reason why programs for constructing rotating rela- 
tivistic neutron star models, like the KEH code [40|, the 
RNS code [Hi, or the BGSM code [H, are not obstructed 
by this nonuniqueness problem is apparently that they all 
utilize an iteration over both the metric and the hydro- 
dynamic equations simultaneously, thereby allowing the 
matter quantities to change during the calculation of the 
metric. 

We want to stress here that these non-convergence is- 
sues in the CFG case are not related to the approximation 
that is made. If one considers this system in the spherical 
(one-dimensional) case, CFC is no longer an approxima- 
tion, but is the choice of the so-called isotropic gauge. 
Even then, the elliptic system P^ - fTT)) no longer con- 
verges to the proper (physical) solution. 



C. The new scheme and its theoretical properties 

Despite the above mentioned convergence problems, 
numerically simulating the physical problem of spherical 
collapse to a black hole in isotropic coordinates has been 
successfully studied by Shapiro and Teukolsky in [ist . 
Because of the spherical symmetry, there exists only one 
independent component of the extrinsic curvature. It 
is then possible to compute directly a conformal extrin- 
sic curvature, t^^K^, from the conserved hydrodynamical 
variables. The elliptic equation for ■!/; then decouples from 
the other elliptic equations by introducing this conformal 
extrinsic curvature and using the conserved hydrodynam- 
ical variables in the source. This source term presents no 
problem for proving local uniqueness, and the equation 
for -0 always converges to the physically correct solution. 
Once the conformal factor, the extrinsic curvature (from 
the conformal factor and the conformal extrinsic curva- 
ture), and the conserved hydrodynamical variables are 
known, the elliptic equation for Nip can be solved and, 
again, the source exhibits no local uniqueness problem. 
This follows from the fact that the extrinsic curvature is 
not expressed in terms of the lapse and the shift. This 
contrasts with the CFC equation ((T6)) where a division 
by N"^ occurs in the last term when the extrinsic curva- 
ture is expressed in terms of its constituents N, ip, and 
In addition, there is no need to use tp' . Finally, the 
elliptic equation for the shift vector can be solved. In 
summary, no problems of instabilities or divergence are 
encountered. 

We now generalize this scheme to the CFC case in three 
dimensions. This involves the use of two different con- 
formal decompositions of the extrinsic curvature: first, 
two different conformal rcscaling and, second, two differ- 
ent decompositions of the traceless part into longitudinal 
and transverse parts. Adopting maximal slicing, K ~ 0, 
a generic conformal decomposition can be written as 

K'^ = V'^-s^^lO^y l^i(ix)'^' + Al^rj^ , (24) 

where C is a free parameter and a a free function, Al^r^^ is 
transverse traceless and L in the conformal Killing opera- 
tor defined by Eq. (fTT|) . We implicitly make use of a flat 
conformal metric, with respect to which Al^rj, is trans- 
verse, although, in principle, it would be more general 
to use the metric 7'-' and the conformal Killing operator 
associated with it, L. But such a decomposition would 
introduce many technical difficulties in our treatment. 
In particular, it is numerically easier to handle tensors 
which arc divergence-free with respect to the flat metric 
in the generalization to FCF. The vector AT*, on which 
L is acting, is therefore called the longitudinal part of 
(yl^(C))u\ xhe first decomposition we use is the one in- 
troduced in Eqs. ([3]) and ([4]) with the choice C = 4 and 
a = 2N. This corresponds to a CTS-like decomposition 
of the traceless part, so that is given by the shift vec- 
tor /?' and Al^rj^ can be expressed in terms of the time 
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derivative of the conformal metric. We denote this trace- 
less part as A'^ := (A^^^^Y^ . In the CFC approximation 
this becomes 

Kv ^ ^-4^y-^ ^ ^{L(3r. (25) 

The second conformal decomposition, 

K'^ = tp-^"A'^, A'^ ^ {LXy^ + Al}^, (26) 

refers to C = — 2 and a = 1. It instead corresponds to 
a conformal transverse traceless (CTT) decomposition of 
the traceless par t of extrinsic curvature introduced by 
Lichnerowicz [4J|. Notice that we have defined A^^ := 
(A(-2))y_ not to be confused with A'^ := (AW)*^'. The 
relation between A'^ and A'^ is given by 

i»i = V^io/r^' =7/,6iy. (27) 

In terms of A'^ , the CFC momentum constraint can be 
written as 

VjA'^ = 8TTi^^°S' = %Tiil^^r^Sj = SttPS*. (28) 

Consistency between the CTT-like decomposition ((26)) 
and the CTS-like one (|25p gcnerically requires a nonva- 
nishing tranverse part Al^j, in Eq. (|26p . However, as it 
is shown in the Appendix, this A^rp is smaller in ampli- 
tude than the nonconformal part /i*-* of the spatial metric 
and A^^ can be approximated on the CFC approximation 
level as 

A'^ w {LXy^ = V'X^ + V^X' - -VkX'T'. (29) 

3 

From Eqs. and an elliptic equation for the vec- 
tor X* can be derived, 

AX' + Iv'VjX^ = SnPS*, (30) 
3 

from which X' can be obtained. With this vector field, 
one can calculate the tensor A'^ via ((29)) . Notice that in 
the case of spherical symmetry, A'"'" = ip^^K" =_d!^K^^ 
is the quantity used by Shapiro and Tcukolsky (isj . 

The elliptic equation for the conformal factor can be 
rewritten in terms of the conserved hydrodynamical vari- 
ables and A'^: 

f-,f- Aim Aij 

AV ^ -27rV'-^^* - V"" g ■ (31) 

This equation can be solved in order to obtain the con- 
formal factor. Once the conformal factor is known, the 
procedure to implicitly recover the primitive variables 
from the conserved ones is possible, the pressure P can 
be computed using the equation of state, and therefore 
S* is at hand. The elliptic equation for Nij} can be re- 
formulated by means of the conserved hydrodynamical 
variables. A'-* , and the conformal factor: 

A(^iV) = 27r7V^-i {E* + 2S*) + 7^"' ^'^"^^'"g^'"'^'' . 

(32) 



From this equation N'>p can then be obtained and, conse- 
quently, so can the lapse function N . Note that, since A'^ 
is already known at this step, no division by N"^ spoils 
the good sign for the maximum principle. 

Using the relation between the two conformal de- 
compositions of the extrinsic curvature, A'^ = tp^A'^ , 
Eq. 1^ can be expressed as {L/SY^ = 2Ni)-^A'K Tak- 
ing the divergence, we arrive at an elliptic equation for 
the shift vector, 

A/3' + ip' {'DjP') = {2N^-^A''^ , (33) 

where the source is completely known. This elliptic equa- 
tion can be solved in order to obtain the shift vector /3* 
consistent with 9t7ij = 0, as required by the CFC ap- 
proximation. 

In this recast of the CFC equations, an extra elliptic 
vectorial equation for the vector field X' is introduced. 
However, now the signs of the exponents of ip and N are 
compatible with the maximum principle for scalar elliptic 
equations, and the problem is linearization stable. While 
this docs not guarantee global uniqueness of the solutions, 
it provides a sufficient result for local uniqueness. This 
strongly relies on the fact that the system decouples in a 
hierarchical way, which we summarize here once more: 

1. With the hydrodynamical conserved quantities at 
hand, solvc Eq. ^ for X\ and thus for A'i . 

2. Solve Eq. (|3ip for t/i, where local uniqueness is now 
guaranteed. Then S* can be calculated consis- 
tently. 

3. Solve Eq. ((5^ for Nip, a linear equation where 
the maximum principle can be applied and unique- 
ness and existence follow with appropriate bound- 
ary conditions. 

4. As the source of Eq. ((55)) is then fully known, solve 
it for I3\ 

Note that this scheme is similar to that used by Shibata 
and Uryu [1^ to compute initial data for black hole - 
neutron star binaries. We will discuss this point further 
in Sec. Ivm 

The new CFC metric equations presented here not only 
allow us to evolve the hydrodynamical equations and re- 
cover the metric variables from the elliptic equations in a 
consistent way (no auxiliary quantity ip' is needed), but 
they also permit to introduce initial perturbations in the 
hydrodynamical variables (strictly speaking, in the con- 
served quantities) in a set of previously calculated initial 
data and directly delivers the correct values for the met- 
ric. It is even possible to perturb only the primitive quan- 
tities, and consistently resolve for the metric by iterating 
until the conformal factor -0, which links the primitive to 
the conserved quantities, converges. We find that such 
an iteration method fails for sufficiently strong gravity if 
the original CFC formulation is used. 
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IV. NUMERICAL RESULTS 

We recapitulate that the original CFC formulation ex- 
hibits serious convergence problems when dealing with 
highly compact configurations such as nascent black 
holes. This weakness of the original formalism is no- 
ticeable in the fact that no simulations of rotational col- 
lapse to a black hole substantially beyond the formation 
of the apparent horizon have been performed so far in the 
CFC. Furthermore, some scenarios which do not involve 
the formation of a black hole are alredy feasible with the 
old formulation only if procedures like using Eq. (j21[) . 
with all associated problems and inconsistencies, are em- 
ployed. An example is the migration of a neutron star 
model from the stable to the unstable branch, which is 
a standard test for relativistic hydrodynamics codes. In 
contrast, the new CFC scheme presented in this work 
solves all problems that prevented performing such sim- 
ulations in the past. In order to show the suitability of 
the new scheme we present the results of numerical simu- 
lations of the migration test and of the rotational collapse 
to a black hole. 



A. Model setup 

The numerical simulations presented here are per- 
formed using the numerical code CoCoNuT [2^, [4^. This 
code solves the evolution of the hydrodynamics equations 
coupled to the elliptic equations for the spacetime met- 
ric in the CFC approximation. Standard high-resolution 
shock-capturing schemes are used in the hydrodynamic 
evolution, while spectral methods are employed to solve 
the metric equations. The code is based on spherical 
polar coordinates, and for the tests presented here we 
assume axisymmetry and symmetry with respect to the 
equatorial plane. Note that the metric equations pre- 
sented in this paper are covariant. Thus the formalism 
can be used for any coordinate basis as well as without 
any symmetry conditions. 

The initial models are general relativistic F = 2 poly- 
tropes in equilibrium with a polytropic constant K ~ 
100. The models are chosen to be situated on the unsta- 
ble branch, i.e. dMADu/dpc < 0, where pc is the central 
rest-mass density. Therefore, any perturbation of the star 
induces either the collapse to a black hole or migration 
to a configuration of the same baryon mass on the stable 
branch. Table U shows the main features of these initial 
models. Models Dl to D4 are uniformly rotating models 
which are identical to those presented in [47j . The model 
labeled SU is a spherical model, while model labeled SS 
is the counterpart model with the same baryon mass but 
it is located on the stable branch. The equilibrium ro- 
tating star models in Dirac gauge (the axisymmetric and 
stationary limit of FCF) used here are described in , 
and are computed using the Lorene [i^ library. We 
map the hydrodynamic and metric quantities to the CFC 
code neglecting the h^^ ~ 10"'^ terms, which arc ncgli- 



TABLE I: Initial models used in the migration test and the 
rotational collapse to a black hole. pc,i is the initial central 
rest-mass density, Qi is the initial angular velocity, rp,i/rc,i 
is the initial ratio of polar to equatorial coordinate radius, 
Madm is the gravitational ADM mass, and J is the total 
angular momentum (which is conserved in CFC during the 
evolution in the axisymmetric case) . Units in which G = c = 
Mq = 1 are used. 



Model 


Pc,i 




rp,i/rci 


re,i 


A/adm 






[10-^] 


[10-^] 










SU 


8.000 





1.00 


4.267 


1.447 





SS 


1.346 





1.00 


7.999 


1.424 





Dl 


3.280 


1.73 


0.95 


5.947 


1.665 


0.207 


D2 


3.189 


2.88 


0.85 


6.336 


1.727 


0.362 


D3 


3.134 


3.55 


0.75 


6.839 


1.796 


0.468 


D4 


3.116 


3.95 


0.65 


7.611 


1.859 


0.542 



gible due to their smallness. Alternatively, we compute 
CFC equilibrium initial models. In this case we find that 
the differences with respect to the FCF models are small 
(~ 0.1%) for representative metric and hydrodynamic 
quantities, initially and during the evolution, and there- 
fore we discuss only the FCF initial models here. 

The hydrodynamic equation are discretized on the fi- 
nite difference grid with x ng grid points. The radial 
grid size is Aro for the innermost cell and increases ge- 
ometrically outwards, while the angular grid is equidis- 
tantly spaced. The metric equations are solved on a spec- 
tral grid consisting of — 1 radial domains distributed 
such as to homogeneously cover the finite difference grid 
and a compactified exterior domain extending to radial 
infinity. On the spectral grid we resolve each radial do- 
main with 33 collocation points. The spherical model 
needs only one angular collocation point, while we use 
17 angular points for the rotating models. 

We track the location of the apparent horizon by means 
of a three-dimensional spectral apparent horizon finder, 
described in detail and tested in [50| • The apparent hori- 
zon location is given by a function 7i(r, 0), which is de- 
composed into a set of spherical harmonics. The coef- 
ficients of Ti. in this basis are computed iteratively, in 
order to satisfy the condition that the expansion in the 
outgoing null direction vanishes at the apparent horizon 
location. 



B. Migration of unstable neutron stars to the 
stable branch 

The first test we consider is the migration of a neutron 
star model in equilibrium from the unstable branch to the 
stable branch, which is a standard but still demanding 
test for general relativistic hydrodynamics codes, as it in- 
volves the dynamic transition between two very compact 
equilibrium states. This test has been performed in the 
past in full general relativistic simulations [5l| . We start 
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the evolution with the nonrotating equiUbrium model la- 
beled SU. Since it belongs to the unstable branch, any 
perturbation from exact equilibrium (which can, for in- 
stance, be caused by discretization errors) leads either to 
a collapse or to an expansion to a new equilibrium con- 
figuration of the same baryon mass on the stable branch. 
The corresponding equilibrium configuration with the 
same baryon mass, model SS, has smaller ADM mass 
than the initial system (see Tabic |T|. Therefore, to pre- 
serve the ADM mass, the final configuration cannot be 
exactly the equilibrium model given by SS. The energy 
difference between models SU and SS should be trans- 
formed into kinetic energy, remaining in the final object 
in the form of pulsations. 

In our case the numerical truncation error is sufficient 
to trigger the migration. Since the final neutron star on 
the stable branch is larger than the initial model (see Ta- 
ble [T|, the outer boundary of the finite difference grid is 
chosen to be 4.5 times the radius of the model SS. We 
perform two simulations on a finite difference grid with 
150 or 300 radial cells and Aro ^ 0.022 or 0.012, respec- 
tively. We use rid = 6 radial domains for the spectral 
grid. Wc evolve the system with either a polytropic or 
an ideal gas equation of state. 

Figure [2] shows the time evolution of the central val- 
ues of the rest-mass density and the lapse. As the star 
expands, pc decreases while Nc grows until the new sta- 
ble equilibrium configuration is reached. In the poly- 
tropic case, there are no physical mechanisms to damp 
the strong pulsations, and the final state resembles a star 
oscillating around the equilibrium configuration until nu- 
merical dissipation finally damps the oscillations. This 
can be seen in the pulsating values of rest-mass density 
and lapse around the value corresponding to the equilib- 
rium model on the stable branch (solid horizontal line in 
Fig.©. 

In the ideal gas case, shock waves are formed at every 
pulsation, and they dissipate kinetic energy into ther- 
mal energy, thereby damping the oscillations. As these 
shocks reach the surface of the star, a small amount of 
mass is expelled from the star and matter is ejected out- 
wards into the surrounding artificial low-density atmo- 
sphere until it leaves the grid across the outer numerical 
boundary. We approximately compute the escape veloc- 
ity as Vc = VWJ w x/ip"^ —• 1, where U is the Newtonian 
potential. This formula is not exact in general relativ- 
ity, but it should by sufficiently accurate near the outer 
numerical boundary where gravity is weaker. We find 
that the shock waves leaving the computational domain 
exceed the escape velocity and therefore the lost mass 
is gravitationally unbounded. Wc also check that these 
results are not affected by changing the resolution or set- 
ting the outer boundary twice as far away. As the oscil- 
lations are damped, the shocks become weaker and the 
mass expelled at each oscillation is smaller. At the end 
of the simulation the star has lost about 10% of its ini- 
tial baryon mass, approaching a state of constant baryon 
mass. As a consequence, the final equilibrium config- 
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FIG. 2: Time evolution of the central rest-mass density pc 
(top panel) and the central lapse A^c (bottom panel) for the 
migration of the unstable neutron star model SU to the stable 
branch, with either a polytropic (solid lines) or an ideal gas 
(dashed lines) equation of state. The dotted horizontal lines 
mark the value of pc and A^c for the equilibrium configuration 
SS from the stable branch with the same baryon mass Mb 
as model SU, while the dash-dotted lines are obtained from 
a series of equilibrium models where mass shedding, like in 
the migration model with an ideal gas equation of state, is 
taken into account. In the inset the baryon mass Mb versus pc 
relation for this model setup is displayed. The models SU (the 
initial model) and SS (the final state for a polytropic equation 
of state) as well as the final state for an ideal gas equation 
of state are marked. The arrows symbolize the respective 
migration paths. 



uration on the stable branch is not the model SS any- 
more but; rather, the corresponding model from the sta- 
ble branch with lower baryon mass and central density. 
In Fig.[2]wc plot the central rest-mass density and lapse 
of a series of equilibrium models on the stable branch 
corresponding to the baryon mass remaining in the com- 
putational domain at each time. It can be seen that these 
values deviate with time from model SS and fit the final 
state in the hydrodynamical evolution of the star. 

As a by-product of this study we draw reader's at- 
tention to the consistency (as it should be) between the 
amplitude and the frequency of the oscillations. The pe- 
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riod of these oscillations is approximately of the order 
of the hydrodynamical characteristic time Tp, which de- 
creases with density like Tp « p^^l"^ . In the polytropic 
case, the maxima of the oscillations in are systemat- 
ically higher than in the ideal gas case. Consequently, 
the characteristic time is shorter than in the ideal gas 
case, as Fig. [5] shows. A second property worth point- 
ing out is that the low numerical viscosity of our code is 
responsible for maintaining a nearly constant amplitude 
of the oscillations (in the polytropic case) during many 
characteristic times. 

Our simulations are consistent with the results from 
a fully relativistic three-dimensional code in [5l| . Sim- 
ilar simulations of this test, with the original, unmod- 
ified CFC scheme, lead to a completely incorrect solu- 
tion with a grossly incorrect ADM mass. When run- 
ning with the new improved CFC scheme, we obtain 
-Wadm = 1.451Afo and initial values for the conformal 
factor and lapse of V'c = 1-561 and = 0.273, respec- 
tively. On the other hand, with the unmodified conven- 
tional CFC scheme, the metric solver already initially 
converges to a solution with Madm = O.647M0 (55%), 
tAc = 1-221 (61%) and ac = 0.532 (63%), where the rela- 
tive differences to the physically correct solution are given 
in parentheses. 

As presented in [H^l the migration test can be suc- 
cessfully simulated using the old CFC scheme, if one re- 
sorts to additionally solving the evolution equation ((2T|) 
for the conformal factor (which would lead to large in- 
consistencies in scenarios with higher compactness but 
still yields acceptable results for the standard migration 
case). Here the superiority of the new, fully consistent 
CFC scheme, which does not depend on such scenario- 
dependend amendments, alredy becomes apparent. 



Collapse of unstable neutron stars to a black 
hole 



As the second test wc present the collapse of a (spher- 
ical or rotating) neutron star model to a black hole. Fol- 
lowing [13] we trigger the collapse to a black hole by 
reducing the polytropic constant K by 2% in the ini- 
tial models Dl to D4. Alternatively, in the spherical SU 
model we increase the rest-mass density by 0.1%, which 
yields a similar dynamic evolution. However, since the 
models are initially in equilibrium, the total collapse time 
depends strongly on the perturbation applied. In these 
cases, the outer boundary of the finite difference grid is 
20% larger than the star radius. For the spherical SU 
model, we perform two simulations using 150 or 300 ra- 
dial cells and Arg ~ 10"'^ or lO^"', respectively, to assess 
the resolution dependence of our simulations. For the ro- 
tating models Dl to D4 the grid is made up of 150 x 20 
and 150 x 40 cells, with the same radial grid spacing as 
in the spherical model. We choose rid = 8 radial domains 
for the spectral grid. As in (47j we use a polytropic equa- 
tion of state in the evolution. 
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FIG. 3: Collapse to a black hole for the spherical model SU, 
and the rotating models Dl and D4. The top panel shows 
the time evolution of the central lapse A^c (thin lines) and 
the central rest-mass density pc relative to the initial value 
Pco (thick lines). The bottom panel shows the time evolution 
of the apparent horizon radius rAH.o in the equatorial plane 
(thin lines) and the rest mass MoutsideAH remaining outside 
the apparent horizon relative to the total rest mass M (thick 
lines). The dashed vertical lines mark the time when the 
apparent horizon first appears. If the axes of the lower panel 
were exchanged, the resulting plot would resemble the typical 
spacetime diagram of a star collapsing to a black hole. 



The top panel of Fig. [2]shows the evolution of the rest- 
mass density and lapse at the center. Since for the max- 
imal slicing condition the singularity cannot be reached 
in a finite time, A^c rapidly approaches zero once the ap- 
parent horizon has formed. In parallel, pc grows, which 
results in a decrease of the numerical time step due to 
the Courant condition applied to the innermost grid cell. 
We terminate the evolution as the central regions of the 
collapsing star inside the apparent horizon become in- 
creasingly badly resolved on the regular grid, and thus 
numerical errors grow. We check in SU model that by 
refining the radial resolution wc are able to follow the 
collapse to even higher densities. Therefore, the only 
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FIG. 4: Isocontours of the rest-mass density for model D4 
after the apparent horizon first appears at t = 129.9. The 
dashed line shows the location of the apparent horizon. 



limitation to perform a stable evolution after the appar- 
ent horizon formation is the numerical resolution used. 
Note, however, that the spatial gauge condition is fixed 
in CFC, and thus we are not able to utilize the common 
method of exploiting the gauge freedom for the radial 
component of the shift vector in order to effectively in- 
crease the central resolution. 

In the bottom panel of Fig. [3] we display the time evo- 
lution of the apparent horizon radius. As expected, the 
apparent horizon appears at a finite radius and already 
encompasses a significant fraction of the total mass of the 
star 70-80%) at that time. Afterwards, its radius 
grows as the surrounding matter falls inside beyond the 
horizon. The fraction of the rest mass remaining outside 
the horizon is also plotted in the figure. In the rotating 
case the apparent horizon is slightly nonspherical. The 
ratio of the polar to the equatorial proper circumferential 
radius of the apparent horizon at the end of the simula- 
tion is i?p/i?e = 0.998-0.978 for models Dl to D4, where 
■^c •= lo'^ y/9^dip/{2'K) and Rp := ^/g^d0/T: 



Since we cannot reasonably determine the location of 
the event horizon, as this would require the evolution 
of spacetime until the black hole has become practi- 
cally stationary, we utilize the apparent horizon radius 
to estimate the mass of the newly formed black hole. 
Following the prescription in [l^l we use the expression 
Mbh = -Rc/2. Note that this formula is only strictly valid 



for a stationary Kerr black hole. In our case, however, 
first of all, some (albeit a small) amount of matter is still 
outside the horizon and the black hole is still dynami- 
cally evolving, and second, the metric of a Kerr black 
hole is not conformally flat [HI]. Still, according to [i^l 
this approximation (excluding the effects of CFC) intro- 
duces an error in the mass estimate of only ~ 2%. For 
the spherical model the estimated value for Mbh at the 
end of the simulation agrees within 0.5% with the ADM 
mass Madm of the initial model, while in the rotating 
models Dl to D3 the error is < 4%. In all these cases the 
above formula overestimates the black hole mass. Be- 
cause of its rapid rotation and the resulting strong cen- 
trifugal forces, in model D4 the collapse deviates signifi- 
cantly from sphericity, leading to a strongly oblate form 
of the density stratification. Consequently, we still find 
a non-negligible amount of matter outside the apparent 
horizon at the end of the simulation (about 12% of the 
total rest mass). Therefore the value for A/bh is 8.2% 
smaller than Madm- In Fig. 3] we present the distri- 
bution of the rest-mass density and the location of the 
apparent horizon at the end of the simulation for this 
particular model. Since the time evolution is limited by 
our chosen, still computationally affordable, grid resolu- 
tion in the central region, we are not able to evolve this 
model to times when a disk forms as in [IJ. Nevertheless, 
all other quantities qualitatively agree with the results in 
that work, although we refrain from performing a more 
detailed comparison due to the respective differences in 
the gauge of the two formulations used in [47| and in this 
study, respectively. 

In the near future we plan to carry out an exhaus- 
tive analysis of the scenario of a collapse to a black hole 
by comparing, on one hand, the CFC formulation with 
FCF (see SecE]), and, on the other hand, the FCF with 
other {free evolution) formulations. The difficulties in- 
duced by the use of different gauges can be overcome 
by using gauge-invariant quantities for comparison and 
analyzing their behavior as a function of proper time. 



V. GENERALIZATION TO THE FULLY 
CONSTRAINED FORMALISM 

The ideas presented in Sec. IIIII can be generalized to 
the FCF approach of the full Einstein equations described 
in Sec. ^TB 

As shown in [^J], the hyperbolic part of FCF can be 
split into a first order system. The reformulation of 
the CFC equations presented in Sec. IIIII relics on the 
rescaled extrinsic curvature A^^ given by Eq. (|27p . Con- 
sequently, we write the FCF hyperbolic part as a first 
order system in {h^^ , A'^), instead of first order system in 
{h^^ , dh'-^ /dt) as in [g^l, arriving at 
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(34) 



(^-^^^^hi^y _ fjkj^tj'^ _ A'^^VkP' - A^^VkP' + "^A^^VkP^ + 2Nil}-^^kiA^^A;>^ 



where 



— i 



(35) 

(36) 
(37) 



The system is closed by adding the equation 

Blip / \ , ■ 2 2 

-gf-T^k [P'wl' + 2Nip-^A'^j = -wlViH^ - f'V„Vif3^ - w'^Vi(3' - ^'^VkVifT + -f^V^Vip' + -w'^Vip\ (38) 

which is derived from applying partial derivatives with respect to t in the definition of w]^ . Moreover, the system 

observes the constraint of Dirac gauge, w^-* = [Eq. ([6])], and for the determinant of the conformal metric, we obtain 
7 = /. The first order system given by Eqs. has the same properties regarding hyperbolicity and existence 

of fluxes as the one in [23| ■ It has the advantage over the second order system for /i*-' proposed in Ref. [13 of getting 
rid of partial derivatives with respect to t of the lapse N, the shift or the conformal factor ip. 
The elliptic part of FCF can be rewritten, using the tensor A^^ , as 



1 



2Trij-^{E* +2S*) 



8ip^ 8 



(NiP), 



(39) 
(40) 
(41) 



r 



The strategy to evolve the two symmetric tensors h^^ 
and A^^ relies on a decomposition of these tensors in lon- 
gitudinal and transverse traceless parts. The longitudi- 
nal parts (divergences with respect to the flat metric) 
are either known a priori or are determined by the el- 
liptic equations. More specifically, the divergence of h^^ 
vanishes according to the Dirac gauge, whereas the di- 
vergence of A''^ is determined by the momentum con- 
straint ((42|) - see below. Consequently, focus is placed 
on the transverse traceless parts of these tensors. The 
latter are described in a pure-spin tensor harmonic de- 
composition, as discussed in a previous article 2J|. In 
particular, each transverse traceless tensor is fully ex- 
pressed in terms of two scalar potentials (named A and B 
in [i^l) that are evolved according to evolution equations 
obtained from the transverse traceless parts of Eqs. 



and ((35l) for h^^ and A'^ , respectively, by applying consis- 
tently the decomposition in [23|. Once the scalar poten- 
tials on the next time slice are determined, the tensors h^^ 
can be reconstructed completely, satisfying the 
divergence-free conditions. This fully fixes h^^ , whereas 
in the case of A*^ the longitudinal part is computed in a 
very similar way to the CFC case, i.e. by determining the 
vector X' from the momentum constraint as described 
hereafter. 

From Eq. (j26|) . the momentum constraint can be writ- 
ten as 



(42) 



which is equivalent to the following elliptic equation for 
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[v'^X^ +V^X*' - '^f^VpXP^ = 

SnfHsn, - f " (v.j^i ~ i^'T- (43) 

This elliptic equation for the vector X* is linear. Since h^^ 
and A'r^rj^ have been calculated previously, we can solve 
the elliptic equation (|43|) to obtain the vector X*. With 
this method, the Dirac gauge and the momentum con- 
straint are guaranteed to be satisfied. Then, A^^ is re- 
constructed from Al^rj, and X^ on the new time slice. 

At this point, since the tensors h'^^ and A*-' are known, 
we can follow exactly the same scheme as in the CFC 
case to solve in a hierarchical way the elliptic equations. 
First the conformal factor is obtained from Eq. ((39)) . then 
the lapse function from Eq. and finally the shift 

vector is acquired from Eq. (|41[) . These equations are 
decoupled in the order mentioned. No sign problems are 
exhibited in the scalar elliptic equation and therefore the 
maximum principle can be applied. A minor concern is 
associated with the sign of the term R in Eq. ([55]) . but 
unique solutions also exist for negative conformal Ricci 
scalars (closely related to R). Note that, contrary to the 
CFC case, here no (additional) approximation has been 
made: it is simply a new scheme to write down FCF, 
where the elliptic part is better behaved from the point 
of view of local uniqueness. Numerical simulations with 
this FCF scheme will be presented in a future publication. 

VI. DISCUSSION 



regimes. We have suggested that these problems are not 
due to the approximative nature of CFC, since FCF in 
the variant of [H, [23| , which is a natural generalization 
of CFC to the nonconformally fiat case, also suffers from 
the same problems. In order to address these issues, first 
focusing on the simpler CFC case, we have considered the 
conformal rcscaling of the traceless part of the extrinsic 
curvature, resulting in the expression for A^^ in Eq. (j27p . 
which is a rcscaling different from the respective ones 
employed in FCF and the CFC approximation, but coin- 
cides with the one in the XCTS approach of [13, El- This 
is motivated by the work of Shapiro and Teukolsky [ist . 
who simulated the collapse of a neutron star model us- 
ing such a reformulation of the CFC metric equations 
(however, restricted to spherical symmetry in their case) 
and apparently did not encounter any of the problems 
described above. Extending their approach to three di- 
mensions, we have decomposed A''^ into longitudinal and 
transverse parts as in the CTT formulation of the con- 
straint equations The divergence (i.e. the longitu- 
dinal part) of this tensor is determined by the momen- 
tum constraints, Eqs. in the CFC case, just as in 
the CTT formulation. In the CFC scheme, we have ne- 
glected the transverse part of this tensor, as the order 
of its error is higher than the one arising from the CFC 
approximation itself. In the nonapproximate FCF case, 
the transverse part of A*-' is determined by an evolution 
equation. Once the conformal extrinsic curvature is ob- 
tained, it can be employed in the Hamiltonian equation to 
calculate the conformal factor tp. The lapse is then fixed 
through the maximal slicing condition, and the resulting 
equation allows the application of a maximum principle 
uniqueness argument. Finally, the shift is found through 
the kinematical relationship defining the extrinsic curva- 
ture, leading to Eq. ([55]) . 



A. Summary 

We have presented an approach to the solution of the 
uniqueness issues appearing in certain constrained for- 
mulations of Einstein equations. We have illustrated the 
problem and its solution through a detailed analytical 
and numerical study of a waveless approximation that 
retains all the involved essential features. 

More specifically, we have reformulated XCTS-like el- 
liptic systems appearing in constrained evolution schemes 
of the Einstein equations, like FCF of [H, [il], as well as 
in the CFC approximation 0, Q . Such systems require 
the simultaneous solution of the constraints, in partic- 
ular, the momentum constraint for the shift, together 
with a maximal slicing condition for the lapse. The re- 
sulting elliptic system presents potential local nonunique- 
ness problems, and numerical implementations have in- 
deed encountered such obstacles. The original CFC for- 
mulation has not been able to cope with these problems, 
as it suffers from convergence of the system to unphys- 
ical solutions or nonconvergence at all in high density 



By performing a variety of tests, we have provided ev- 
idence that the problem of convergence to an unphysi- 
cal solution of the metric equations (or even complete 
nonconvergence) in the original formulation of the CFC 
scheme is fully cured by our new reformulation. Not 
only can numerical results in the original CFC scheme (in 
the, at most, moderately gravitationally compact regime 
where that system still yields physically correct solutions) 
be reproduced by the new formulation but, more im- 
portantly, the new numerical results presented here ex- 
hibit the proper numerical and physical behavior even for 
highly compact configurations. For the first time, it has 
been possible to successfully perform both the migration 
test and the collapse of a neutron star to a black hole in 
the CFC case in a consistent way. Our new formulation 
thus facilitates simulations in the high density regime of 
those scenarios where the CFC is still a reasonably fair 
approximation, that is, for systems which are not too far 
from sphericity, like stellar gravitational collapse. 
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B. Comparison with previous works 

As compared to the original CFC formulation by Iscn- 
berg Q and Mathews and Wilson [5], the scheme pre- 
sented here is augmented by an additional vector elliptic 
equation for X', while the elliptic character of the system 
of metric equations is preserved. The new scheme refor- 
mulates the CFC approximation in a CTT shape (one 
scalar and one vector elliptic equation), and then solves 
for the lapse and the shift (one additional scalar and one 
vector elliptic equation). In contrast, the original CFC 
scheme employed an (X)CTS approach where, together 
with two scalar elliptic equations, only one vector elliptic 
equation was present. In contrast to the original scheme, 
the elliptic system in the new formulation not only cor- 
rects the problem of local uniqueness in the scalar ellip- 
tic equations, but also introduces a hierarchical structure 
that decouples the system in one direction. 

In the context of the conformally flat approximation, 
the same "augmented CFC" scheme as that discussed 
here has been introduced already by Saijo [1] to compute 
gravitational collapse of differentially rotating supermas- 
sive stars. However, in this work the inconsistency be- 
tween Eq. (pSj) and Eq. i.e. setting to zero the trans- 
verse traceless part of A^^ , has not been pointed out. On 
the contrary, we have analyzed this inconsistency in de- 
tail (cf. the Appendix) and have shown that it leads to an 
error of the same order as that of the CFC approximation. 
In addition, we have shown here that the introduction of 
the vector potential is the key ingredient for solving 
the nonuniqueness issue. 

The same scheme, but without the conformal rescal- 
ing of the matter quantities, has also been used recently 
by Shibata and Uryii 45[ in the context of computing 
initial data. As in the inconsistency resulting from 
setting to zero the transverse traceless part of A^^ and 
the uniqueness issue are not discussed in their work. We 
emphasize that the these studies 0, |4^ do not discuss 
the extension of the new scheme to the nonconformally 
flat case, as done here. 

Let us also mention that the augmented CFC scheme 
presented here can be regarded as a hybrid mixture of 
some of the waveless approximation theories (WAT) pro- 
posed by Isenberg [1]. In fact, the CFC approximation 
using the two choices 7y ~ fij and dt'^ij — [as employed 
in Eq. p3|) ] corresponds to version WAT-I. On the other 
hand, the approximation Al^rj, = used in Eq. is 
in the spirit of the vanishing transverse traceless part of 
the extrinsic curvature in the (coupled) version, WAT-II 
(although WAT-II refers to the physical extrinsic curva- 
ture, whereas here we have dealt with the conformal one). 
As mentioned above, both assumptions are consistent at 
the considered level of approximation, as shown in the 
Appendix. 

Regarding the complete constrained evolution of the 
Einstein equations, we have generalized the ideas pre- 
sented here for the CFC case to the elliptic part of FCF. 



In previous studies [23|, |2J] , the hyperbolic part of Ein- 
stein equations resulted in a wave-type equation for the 
tensor ft,'-' , representing the deviation of the three-metric 
from conformal flatness. With the introduction of A^^ 
we have recovered here a first-order evolution system, 
analogous to the standard Hamiltonian 3 4-1 system, in 
which we have, however, retained only the divergence- 
free terms. Thus, for both h^^ and A^^ , the transverse 
(divergence-free) parts are evolved by this system, while 
the longitudinal parts are fixed either by the gauge (for 
/i*-' ), or by the momentum constraint (for A^^). Numeri- 
cal results for this case will be presented in future studies. 

We finally comment on the recent work by Rinne [3 , 
where uniqueness problems appearing in certain con- 
strained and partially constrained schemes for vacuum 
axisymmetric Einstein equations [20l . [s^ are addressed. 
As in the present case, uniqueness issues related to the 
Hamiltonian constraint equation are solved by adopting 
an appropriate rescaling the extrinsic curvature. On the 
other hand, problems associated with the slicing condi- 
tion are tracked to the substitution in that equation of 
the extrinsic curvature by its kincmatical expression in 
terms of the (shift and the) lapse. The latter spoils the 
uniqueness properties by reversing the sign of the rele- 
vant term in the slicing equation. This problem is solved 
by enlarging the elliptic system with an additional vec- 
tor so as to reexpress the relevant components of the 
extrinsic curvature without resorting to the lapse. The 
resulting elliptic system also presents a hierarchical struc- 
ture. Although the spirit of such approach is close to 
the one here presented, the specific manner of introduc- 
ing the additional vector variable in [l^ critically relies 
on the two dimensionality of the axisymmetric problem 
(specifically, on a choice of a particular gauge and on 
the fact that vectors and rank-two traceless symmetric 
tensors have the same number components in two di- 
mensions, a property lost in three dimensions). On the 
contrary, the introduction of the vector X' through the 
CTT decomposition ([^ is properly devised to work in 
three dimensions. Relevant related discussions in the 
three-dimensional context can be found in Sec. 3.4 of [l^ 
(where the relation between nonuniqueness problems in 
XCTS and axisymmetric constrained evolution schemes 
is discussed) and in the three-dimensional constrained 
evolution scheme presented by Moncrief et al. in [25| . 
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APPENDIX A: CONSISTENCY OF THE 
APPROXIMATION 

In the derivation of the new formalism we make use 
of the fact that {LXy^ w A'^ in CFG. We show next 
that this assumption is completely consistent at the ac- 
curacy level of the CFG approximation. In the first place, 
we need to estimate the error of the CFG approximation 
itself. By definition, the CFG three-metric deviates lin- 
early with /i*-' from the (exact) FGF case. It can be easily 
shown from the FGF equations (|5^ - (PT|) that the metric 
quantities behave as 

^ = ^CFC+OW, (Al) 
N = NcFC+Oih), (A2) 
[3^ =PhFC+0{h). (A3) 

Therefore h^^ can be used as an estimator for the error 
of the CFG approximation. 

Two limits in which CFG is exact will be considered. 
First, in spherical symmetry the CFG metric system is an 
exact reformulation of the Einstein equations since /i*-' = 
in the FGF metric. If the system is close to spherical 
symmetry (i.e. spheroidal), and if we are able to define 
a quasispherical surface of the system (e.g., the surface 
of a star or the apparent horizon of a black hole), then 
the equatorial and polar circumferential proper radius. 
Re and Rp, can be computed, and we can define the 
ellipticity of the system as 

1 - Rl/Rl (A4) 

Close to sphericity, scales linearly with h^^ , and we 
can ensure that the error of CFG is h^^ ^ C'(e^). The 
second limit to consider is if a post-Newtonian expansion 
of the gravitational sources is possible, i.e. if the post- 
Newtonian parameter ma,x{v^ /c^,GM/Lc^) < 1, where 
17, Af , and L are the typical velocity, mass, and length 
of the system, respectively. In this case the CFG metric 
behaves like the first post-Newtonian approximation (ssl . 




V' = ^CFC + O (l/c^) , (A5) 
N = NcFC + O (l/c^) , (A6) 
cP' = cPhFC + 0{l/c*). (A7) 

Note that, for clarity, we explicitly retain powers of the 
speed of light c as factors in the equations throughout 



this appendix. In the case that both limits arc valid, i.e. 
close to sphericity and in the post-Newtonian expansion, 
the nonconformally-flat part of the three-metric behaves 
like h"^^ ~ 0{e'^/c'^). The next step is to compute the 
behavior of the CFG metric if we assume {LXy^ « A^^ , 
considering the two limiting cases introduced above. 

In the spherically symmetric case the relation 
{LXy^ ^ A'i is trivially fulfilled. Therefore the behav- 
ior for a quasi-spherical configuration is also /i*^ ^ C'(e^) 
even if Al^r^^ = is assumed. This limit in the approxi- 
mation is very important, since it is independent of the 
strength of the gravitational field. For example, it allows 
us to evolve black holes, with the only condition that h^^ 
should be small, i.e. close to the sphericity. 

To check the approximation in the post-Newtonian 
limit, we need to compare /3cfc ^^'^ This can be 
done by means of the post-Newtonian expansion of the 
sources of Eqs. pT)) and ((30)) . respectively, 

A/?^FC + I'D'T^jPhFC = IStt^*' + O (1/c^) , (A8) 

AX' -I- -V'VjX' = SttS" + O (1/c^) . (A9) 
3 

From the comparison of Eqs. (|A8[) and (|A9p we obtain 
that 

c^^ = c^X' + 0{l/c^). (AlO) 
Thus A^^ can be computed in terms of X' as 

c'A^' = ^^c\LPcFcr = c\LXr + O {II c^) , 
^jVcfc 

(All) 

where we make use of V'cfc/^CFc = 1 -I- C(l/c^). The 
effect of using {LXy^ instead of A^^ in the calculation of 
the CFG metric can be seen in the expressions 

iJCFC ^ A-^S(^^)(NcFC,i'CFC,A'^) 

= A7i5(^)(iVcFC, V'CFC, (LXy^) 

+0 (l/c«) , (A12) 

A^CFC = l/'cFcAs"^'5(W'/0(^CFC,V'CFC,^*^') 

= ^cFcAs"''5(w^)(A^cFC,V^CFC, (LXy^) 
+0 (l/c8) , (A13) 

c/9cFc = cA-i5(^)(7VcFC,V'CFc,i'^') 

- cA;l5(^)(iVcFC>CFC,(W') 

+0 (l/c^) . (A14) 

where iS(^), 5(jvi/;) and iS(^) are the sources of Eqs. (pij) - 
((33)) . and A^^ and A^^ are just the inverse operators ap- 
pearing in the right-hand-side of these equations (for the 
scalars ijj and Nip, and for the vector respectively). 
When comparing Eqs. (|A12p - (|A14p with Eqs. ((X5|) - ((X7)) . 
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FIG. 5: Consistency of the approximation for rotating neu- 
tron star models. In the top panel max |j4^rp/y4*-' | for FCF 
(solid line) and CFC (dashed line) as well as the maximum 
deviation from conformal flatness max|/i*^ | for FCF (dash- 
dotted line) are plotted against the ellipticity e. The bottom 
panel shows the absolute difference |A'^c,cfc — A'^d in the cen- 
tral value of the lapse between CFC and FCF (solid line) 
and the absolute difference |A^c,cfc — N^^qyc I between regu- 
lar CFC and CFC neglecting A^^ in Eq. ^ (dashed hue). 
The Kepler limit is marked by vertical dotted lines, while 
the slanted dotted lines represent the order of accuracy with 
respect to powers of e. 

it becomes obvious that in all cases the error introduced 
by making the approximation {LXy^ « A^^ is smaller 
than the error of the CFC approximation itself. 

As an illustration of the above properties, we study the 
influence of the A!^j, term in Eq. (j29p when computing 
rotating neutron star models with a polytropic F = 2 
equation of state. This model setup contains the initial 



models used in Sec. IIVI They assume axial symmetry 
and stationarity, in combination with rigid rotation. We 
build a sequence of rotating polytropes with increasing 
rotation frequencies, while keeping the central enthalpy 
fixed, which produces models of increasing masses from 
AI = 1.33 Mq (no rotation), to M = 1.57Mq (the Ke- 
pler limit; see below). For all these models, we use three 
gravitational field schemes: the exact Einstein equations 
using the stationary ansatz in FCF, and the two approx- 
imate ones, regular CFC and CFC, neglecting the term 
A^j^rp in Eq. ([^5]) . The results are displayed on a logarith- 
mic scale in Fig. [J] In the top panel we show the maximal 
amplitudes of Al^j, (relatively to A^^) in both FCF and 
regular CFC, as functions of the ellipticity e defined in 
Eq. (jA4p . This quantity is physically and numerically 
limited by the minimal rotational period at the so-called 
mass-shedding limit (or Kepler limit), when centrifugal 
forces exactly balance gravitational and pressure forces 
at the star's equator. In the FCF case we plot the maxi- 
mal amplitude of h^^ . This quantity is dimensionless and 
represents the deviation of the three-metric from con- 
formal flatness, which can be interpreted as the relative 
error one makes in the metric when using CFC instead of 
FCF. Note that this error in computing A''^ by discard- 
ing the A^rp term in the CFC approximation is roughly 
of the same magnitude as the error on the metric in the 
CFC case. All these quantities decrease like 0{e^) as ex- 
pected, except for stars rotating close to the Kepler limit. 
Indeed, the development in powers e is equivalent to a 
slow- rotation approximation (see, e.g., [s^ ) by perturb- 
ing spherically symmetric conflgurations, and, when com- 
paring these slow-rotation results to numerical "exact" 
ones for rigidly rotating stars (see, e.g., [Ill in the two- 
fluids case), one sees that they usually agree extremely 
well, excepted very close to the Kepler limit, where this 
"perturbed spherical symmetry" approach is no longer 
valid. Finally, because A^^ appears as a quadratic source 
term in the Poisson-like equations (fT5|) and , the over- 
all errors on the lapse or the conformal factor ijj are 
even smaller, as shown in the bottom panel of Fig. [5l In 
the case of the central value iVc of the lapse, the error due 
to the CFC approximation is maximal at the Kepler limit 
and < 10~^ for the studied sequence. The error which 
is then due to neglecting Al^j, within the CFC scheme 
amounts to < 10~^ and decreases faster than the error 
due to the CFC approximation, namely, as 0{e'^), again 
except near the Kepler limit. Our tests thus show that for 
stationary rotating neutron star models this additional 
approximation induces an error which falls within the 
overall CFC approximation. 
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